Steady Periodic Shear Flow is Stable in Two Space Dimensions . 

Nonequilibrium Molecular Dynamics versus Navier-Stokes-Fourier 

Stability Theory — A Comment on two Arxiv Contributions . 

Q ■ Wm. G. Hoover and Carol G. Hoover 

cni : 

5_i . Ruby Valley Research Institute 

^ I Highway Contract 60, Box 601 

^ ■ Ruby Valley, Nevada 89833 

^ : 

(Dated: February 25, 2013) 

S^ ■ Abstract 

<D • 

^ . Dufty, Lee, Lutsko, Montanero, and Santos have carried out stability analyses of steady sta- 
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C^ , tionary shear flows. Their approach is based on the compressible and heat conducting Navier- 
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) ' ' Stokes-Fourier model. It predicts the unstable exponential growth of long-wavelength transverse 
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S ' perturbations for both two- and three-dimensional fluids. We point out that the patently-stable 

'^ , two-dimensional periodic shear flows studied earlier by Petravic, Posch, and ourselves contradict 

O I these predicted instabilities. The stable steady-state shear flows are based on nonequilibrium molec- 

ular dynamics with simple thermostats maintaining nonequilibrium stationary states in two space 
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K* . dimensions. The failure of the stability analyses remains unexplained. 

o: 

en , PACS numbers: 05.20.Ji, 47.ll.Mn, 83.50.Ax 

(T^ I Keywords: Molecular Dynamics, Shear Instability, Navier-Stokes-Fourier 



X 



I. INTRODUCTION 

It recently came to our attention that the stability of high-Reynolds'-number atomistic 
shear flowsi"- apparently contradicts the instability predicted by a perturbation analysis^*^ of 
the Navier- Stokes- Fourier equations for steady compressible heat-conducting simple shear. 
Here we develop an alogorithm for solving the Navier- Stokes- Fourier problem directly. We 
also apply a two-dimensional version of the perturbation analysis to soft disks. This work 
confirms the contradictory nature of the microscopic and macroscopic approaches. The 
failure of the perturbation analysis remains unexplained. 

In Section II we describe the atomistic and continuum approaches to simulating simple 
shear. We apply the continuum algorithnt^i^ to a Navier- Stokes- Fourier gas-phase model 
of simple shear. In Section III we formulate a dense-fluid constitutive model^'^'ii for the 
continuum description of the soft disks used in atomistic simulations. In Section IV we apply 
several versions of the continuum perturbation theory to this soft-disk constitutive model. 
The results of that theory are contradicted by the existing molecular dynamics results. 
Section V sums up the current state of our knowledge, emphasizing the open problem which 
remains unsolved. 

II. STEADY SHEAR FLOWS IN TWO SPACE DIMENSIONS 

Steady homogeneous shear flows with periodic boundaries would seem to be the simplest 
conceivable nonequilibrium steady states. Figure 1 shows a particle snapshot taken from a 
typical molecular dynamics simulation. In that Figure the central cell, containing A^ = 100 
two-dimensional particles, is repeated periodically in space. Periodic images to the right 
and to the left are generated by adding or subtracting multiples of the cell side length L to 
all the particle coordinate values from the central cell. 100-particle images just above the 
central cell move to the right with an additional velocity +eL relative to the central row of 
cells. Images just below that row likewise move to the left with an additional velocity —eL. 
These periodic boundaries impose a shear strain rate and a corresponding stress on all cells. 
The strain rate e defined and imposed by these periodic boundary conditions is the overall 
macroscopic velocity gradient : 

e = ( du^/dy ) . 



N = 100, Periodic Shear Flow 




Figure 1: Snapshot from a molecular dynamics simulation of steady periodic shear in two space 
dimensions. The three periodic images of the central L x L cell shown in the top row move to the 
right with velocity +eL relative to the central 100-particle cell. The three images in the bottom 
row move to the left with velocity — eL . The regular checkerboard arrangement of the cells recurs 
at integral multiples of the shear period, r = (1/e) . 



Such early 100-particle studies from the 1970s were, by the 1990s, augmented by larger- 
scale simulations with as many as 514^ = 264, 196 particles^"-. These steady-state flows, in 
which the heat generated by the shearing motion is extracted by "thermostat" (or ergostat) 
forces, established that the measured shear viscosity, at fixed strain rate, is insensitive to 
system size. 

Consider a prototypical dense soft-disk fluid, with N = L"^ and a short-ranged repulsive 
pair interaction 0(r < 1) = 100(1 — r^)^. For A^ greater than 64 the shear viscosity has only 
a small "size dependence". For a fixed strain rate the finite- A^ viscosity ri{e,N) approaches 
the large-system limit from below. The iV-dependent deviations from the limiting viscosity, 



ri{e,N — > cxd) , are of order J{l/N) . The deviation is about two percent for A^ = 64 and is 
negligible, relative to statistical fluctuations, for systems with A^ of order 10, 000^"^. 

One might well expect, by analogy with macroscopic three-dimensional hydrodynamic 
flows^^, that a turbulent instability would appear at a Reynolds' Number of the order of a 
thousand, R = eL? /{tj/ p) ~ 1000 . In fact, detailed investigations, for thermostated two- 



dimensional flows with Reynolds' Numbers as large as 50,000 , showed no instability what- 
ever. These molecular dynamics simulations showed instead that steady two-dimensional 
dense-fluid shear flows are stable. These results were announced in 1995 in two papers by 
author Bill and Harald Posch, "Large-System Hydrodynamic Limit" and "Shear Viscosity 
via Global Control of Spatiotemporal Chaos in Two-Dimensional Isoenergetic Dense Fluids" . 

Soon afterward two papers suggesting a long-wavelength exponential instability of all 
such flows appeared in the Condensed Matter arXiv: "Stability of Uniform Shear Flow" 
and "Long Wavelength Instability for Uniform Shear Flow" , by Jim Dufty, Mirim Lee, Jim 
Lutsko, Jose Montanero, and Andres Santos^'^. These papers considered linear perturbations 
of the density, velocity, and energy flelds from the standpoint of the Navier- Stokes- Fourier 
equations, specialized to small strain rates e and to long wavelengths (or small k vector) 
normal to the flow direction (A = 2Ti/k) . Their analysis showed that for any flxed strain 
rate exponential growth invariably results for perturbations of sufficiently long wavelength. 

In early 2012 we chanced upon these contradictory results [ stable molecular dynamics 
and unstable linear perturbation theory, both for the same problem ] . To address this 
puzzle here we first investigated solutions of gas-phase Navier-Stokes-Fourier shear fiows 
with moving periodic images of a square Eulerian continuum grid, using the same periodic- 
shear boundary conditions illustrated in Figure 1 . 

A. Continuum Algorithm for Steady-State Simple Shear 

The finite-difference numerical work-i^i^ is most naturally carried out for grid space and 
time differences Ax and At chosen near the Courant condition limit. At < Ax/c , where c 
is the fiuid's sound speed . The equations to be solved are minor modifications of the usual 
conservation laws for mass, momentum, and energy : 

p = -pV -u] pu = -V ■ P ; pe = -P -.Vu-V -Q . 

The boundaries, where the square middle-row cells of Figure 1 contact the upper and lower 
rows of cells, need to be made consistent with steady shear. Here p is the mass density, 
u = {ux.Uy) is the velocity, P is the pressure tensor, e the energy per unit mass, and Q 
the heat fiux vector. The superior dots indicate "comoving" "Lagrangian" time derivatives 
following the motion. An efficient algorithm solving this problem can be based on those 



developed for the Rayleigh-Benard continuum convection problems described in References 
5, 8, and 9. 

In the continuum periodic-shear algorithm we added a "global" ergostat to the energy 
equation to constrain the total internal energy. We also added two global "thermostats" to 
the continuum equations of motion. These thermostats constrained the summed-up squares 
of the X and y velocity fluctuations to small constant values, relative to the mean simple- 
shear motion Ux{y) = ey . We also included a short-ranged "artiflcial viscosity" to avoid 
an "even-odd" instability. The one-dimensional analog of this instability would lead to 
alternating signs for velocities evaluated at neighboring nodes of a periodic computation 
mesh with A^ zones deflned by A^ nodes. At the end of each timestep the ith nodal velocity 
is replaced by a weighted average (0 < a < 1) including the zone velocities of its neighboring 
zones : 

u-^au^ + {l-a){ul + uU)/2. 

The two-dimensional periodic shear problem is implemented with densities evaluated 
in A^ computational square Aa; x Aa; zones with velocities and energies evaluated at the 
nodes deflning the four corners of the zones. Fourth-order Runge-Kutta integration of the 
4 X A^ = 4L^ ordinary flnite-difference equations results in a natural strain rate , 

e = 2Aa;/(LAt) . 

The factor 2 reflects the two increments of time occuring in a single Runge-Kutta time 
step. This strain rate is of the same order as the k vector describing the longest-wavelength 
perturbation : 

k = {271/1) ~ 7re(At/Aa;) ~ vre/c , 

and is also fairly close to the critical wavelength of the instability analyses. 

B. Typical Numerical Results for Continuum Shear of a Gas 

A typical 512 x 512 zone simulation, with 

[ Ax = 1 ^ L = 512 ; At = 0.25 -^ e = 2Aa;/(LAt) = (1/64) ] ; 
[ (PV/NkT) = (E/NkT) = l^c=V2; u = {r]/p) =p = r] = K=l; r]v = 0], 
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Figure 2: Time history of the horizontal kinetic energy, which includes the single-cell contribution 
of the shear to the kinetic energy (below) and the much smaller vertical kinetic energy (above) in 
a square 32 x 32 Eulerian unit cell. The vertical kinetic energy scale has been multiplied by 10^ . 
The strain rate is e = (2Ax/At)/L = (2/(l/2))/32 . A portion of the time history corresponding 
to strains from 900 to 1000 is shown. The transport coefficients, density, and internal energy, 
{ ri,K, p,e }, are all equal to unity. The Reynolds Number is i? = 4 x 32 = 128 . 

corresponds to a shear flow with a Reynolds' number 

R = L^e/z/ = 512 X (512/64)/z/ = 4096 , 

and is perfectly stable, for a run time of 256 x 50 steps of 0.25 each, corresponding to a total 
shear of 50 . 

We carried out a variety of simulations. Time-dependent data from one of them are 
shown in Figures 2 and 3 . The kinetic energy history, after an irregular transient, is typical, 
with nearly periodic oscillations of the horizontal and vertical velocity deviations from simple 
shear flow [ these deviations are rescaled at the conclusion of each timestep ] . The deviations 
are analogous to thermal fluctuations. Figure 3 shows these same fluctuations as arrows. 
These patterns are not stationary, but vary with time. Because the time periodicity of 
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Figure 3: The velocity, relative to the mean flow, x = ey for the 32 x 32 grid described in Figure 2. 
Although this velocity perturbation is not stationary - it changes with time - this view is typical. 

the shear, with period r = (1/e), is built into the boundary conditions an inhomogeneous 
perturbation flow is to be expected. 

We explored systems with L x L computational zones, with L = 8, 16, 32, 64, 128, 256, 
and 512 to strains of at least 100. With Ax = 1 and At = (1/2) , so that the Reynolds 
number is R = 4L , all members of this series were stable, giving results similar to those 
in Figures 2 and 3 . Somewhat larger or smaller values of At typically led to instability, 
making it difficult to validate the stability work of References 6 and 7 . 

These exploratory calculations could be extended in many arbitrary ways. We chose 
instead to follow and extend the linear perturbation-theory hard-sphere work of Dufty et 
alii^'''' directly to a numerical model based on our earlier soft-disk shear-flow model. Our 
results, described in what follows, suggest that the linear perturbation analysis yields flawed 
predictions, contradicted by the well-established stability of thermostated simple shear in 
two space dimensions'-^^. Despite our efforts to extend the perturbation analyses, described 
below, a convincing theoretical demonstration of this finding remains a challenging research 
goal. It seems likely that the flaw in the perturbation treatment is an inadequate treatment 



of steady-state thermostating. Our checks of the other aspects of the perturbation theory 
have so far revealed no errors. 

III. CONSTITUTIVE MODEL FOR SOFT DISKS 

The short-ranged repulsive potential 0(r < 1) = 100(1 — r^)^ has been carefully inves- 
tigated near the reference state, [ p = 1 ; e = 1 ] , where e is the energy per particle and 
each particle has unit mass^. The equilibrium temperature and pressure at that state are 
approximately T = 0.69 and P = 3.90 . For our numerical stability studies we adopt the 
equation of state from Reference 3 : 

P/p = 5 + 8(p - 1) + 2.5(e - 1.443) + 9(p - 1)^ + 2(p - l)(e - 1.443) ; 

T = 1 - (p - 1) + 0.7(e - 1.443) - 0.8(p - 1)^ - 0.5(p - l)(e - 1.443) ; 

e - 1.443 = 1.5(p - 1) + 1.5(r - 1) + 2.4(p - 1)^ + 1.2(p - l)(r - 1) . 

We estimate the transport coefficients (and evaluate their density and temperature deriva- 
tives numerically) rj, rjv, k from Enskog's theory^^'^^. 

Enskog suggested that transport properties for soft potentials be estimated from those 
of a corresponding hard-potential model. The equivalence is based on setting the soft 
model's "thermal pressure", T{dP/dT)y, equal to the hard-potential pressure. In our two- 
dimensional implementation the hard-disk model corresponding to the soft-disk potential 
= 100(1 — r^)^ follows from the relation : 



{V/Nk){dP/dTUt = (PV/NkT) 



hard 



For the hard-disk compressibility factor we use the ten-term virial series from Clisby and 

McCoy's workia : 

[ {dP/dT)/p ]soft = iP/pT)t.rd = 1 + ibp)x = l + bp + 0.782{bpf + 0.532(6p)3+ 

0.334(6p)^ + 0.199(6p)^ + 0.115(6p)*^ + 0.065(6p)^ + 0.036(6p)^ + 0.020(6p)^ + . . . . 

Here x is the ratio of the actual hard-disk collision rate from the virial series to the low- 
density kinetic-theory collision rate calculated with the two-term series. The computed value 



of the hard-disk second virial coefficient b = 7ro"^/2 gives the hard-disk diameter a = J2b/7i 
needed for the hard-disk transport coefficients. 
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David Gass kindly worked out formulae for the Navier-Stokes-Fourier transport coeffi- 
cients for hard disks in 1970^^ : 



7] = 7]obp[ (bpx)-' + 1.0 + 0.8729{bpx) ]; Vv = l.2AQr^ohp{hpx) ; Vo = JmkT/7i/2a 



K = Kobpl {hpxY^ + 1.5 + 0.8718(6px) ] ; '^o = 2^k^T/mTx/a . 

These expressions provide estimates for the soft-disk transport coefficients needed to apply 
the perturbation theory described in the following Section. 

IV. DUFTY-LEE-LUTSKO-MONTANERO-SANTOS' LINEAR THEORY 

Consider the reference steady shear flow state for soft disks at unit density and energy : 

{p=l, U:c = ey , Uy = 0, e= 1 } . 

There are extensive molecular dynamics simulation data available for strain rates e of 0.10, 
0.25, and 0.50 . We apply linear perturbation theory to this problem for a perturbation with 
components 

{ 6p X e'^'^+*"* , 6u^ X 6'^=^+*^* , 6uy x 6'^^+'"'' , 6e x e*'=^+''^* } , 

where the unperturbed state is steady shear at a mass density and internal energy per 
particle of unity : 

{ 5p = p-1 , 6ux = u^- iy , 6uy = Uy , 6e = e- I } . 

In addition to the conservation laws we use the Navier-Stokes-Fourier constitutive relations 

P = [ Peq - AV ■ M ]/ - r/[ Vm + Vm* ] ; Q = -kVT , 

with the Gass-Enskog transport coefficients. As usual rj is the shear viscosity. Here A is the 
difference between the bulk and shear viscosities , X = rjv — ri , and k is the heat conductivity. 
To a good approximation A ~ 0, so that the bulk and shear viscosities are roughly equal, 
and about four times less than the heat conductivity : 



A. Matrix Solution of the Linear Equations 

If we represent the hydro dynamic state perturbations by the four- dimensional vector 

5 then the hnearized hydrodynamic equations take the form of a 4 x 4 matrix equation, 

6 = M ■ 6. Let us begin by estimating the values of the 16 matrix elements. If we (1) ignore 
the nonlinear convective contributions to the equations of motion and (2) assume a solution 
of the form exp[ iky + iut ] we get a set of four linear equations for the evolution, in space 
and time, of perturbations in mass, momentum, and energy : 

[ mass conservation ] : p = — pV ■ u — > 

{d6p/dt) = iu6p = —p{d6uy/dy) = —ikpSuy ; 

[ X momentum conservation ] : iix = —e6uy — {1/ p){dPxy/dy) — )■ 

iu6u^ = -e6uy - {l/p){d/dy)[ -r]{d/dy){ey + 6u^) ] ~ 

-e5uy + ik{e/p)[ {dr]/dp)e5p + {dr]/de)p5e ] - {ri/p)k^5ux . 
[ y momentum conservation ] : Uy = —{l/p){dPyy/dy) — > 

iuSuy = -{l/p){d/dy)[ Peg - (v + r]v){duy/dy) ] ~ 

-{ik/p)[ {dPeg/dp)e6p + {dPeg/de)p6e ] - [r] + r]v)k'^6uy . 
[ energy conservation ] : e = (— l/p)[ P : Vu + V ■ Q ] = 

-(1/P)[ P,y{dujdy) + Pyy{dUy/dy) + (dQy/dy) ] ^ 

iuj5e ~ e^[ {dv/dp)e5p + {dv/de)p5e ] + 2ikvt5ux 

-ik{Peg/p)6uy - P{k/p)6T . 

Note that the kinematic viscosity, u = {fj/ p) has been introduced in the energy conservation 
expression. For computation of the matrix elements, and their eigenvalues, it is necessary 
to express V^T in terms of 5p and 5e : 

5T = -0.786p + 0.76e -^ V^T ~ k'^[ 0.786p - 0.76e ] . 

A numerical evaluation of the eigenvalues for a dense grid of A;- vectors and strain rates gives 
two unstable ( positive real part, leading to exponential growth ) eigenvalues in the shaded 
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Figure 4: Portion of /c-vector versus strain rate space for two-dimensional soft disks. Flows in a 
lower triangular region [ shading shows that region for the Sllod Matrix ] are unstable, leading to 
longtime exponential growth of perturbations, according to the theory outlined in References 6 and 
7. Stable atomistic simulations for soft disks have been carried out for the conditions corresponding 
to the three filled circles (as well as for many other situations). The two approaches to stability 
give contradictory results. The straightforward nature of the molecular dynamics argues that the 
perturbation-theory results are invalid. 

lower-triangular region shown in Figure 3 ( which closely resembles the behavior found by 
Dufty et alii ). Notice that the stable large-system molecular dynamics data from References 
3 and 4, shown also in the Figure, lie squarely inside the "unstable region" . 

The contradiction between perturbation theory and molecular dynamics is robust. We 
have investigated several modifications of the basic theory and found qualitatively similar 
results. Replacing the "Sllod" matrix element with its "Doll's- Tensor" analog, 

{ M^^^uy = -e ; M,^,„, = o } ^-^ { M„,,„^ = ; M,^,„, = -e } , 

leaves the qualitative picture unchanged. Reducing the matrix to 3 x 3 by ignoring the 
possibility of energy fluctuations still yields an unrealistic triangular "unstable region" . We 
have no explanation for the failure of perturbation theory to deal with steady simple shear. 

V. CONCLUSION 



Although the evidence from molecular dynamics supports the large-system stability of 
thermostated steady shear, linear stability analysis, with a fixed global thermostat, predicts 
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instability. Where might the explanation he? We don't know. Explanations might involve 
either of two flow features treated poorly by the theory: [1] The "steady-state" shear flow 
is not really stationary. See again Figure 2. At times which are integral multiples of the 
inverse shear rate, t = n{l/e), the periodic arrangement of cells is a perfect rectangular 
checkerboard; [2] Evidently the molecular dynamics thermostats provide additional stability 
to the flow. 

A variety of molecular dynamics algorithms, with thermostats or with ergostats and with 
Doll's or Shod equations of motion, all give very similar results, all of them identical for 
the linear response proportional to e . It may be that some more faithful representation 
of the thermostat forces in the perturbation theory could lead to results consistent with 
the molecular dynamics. We would be most grateful for an explanation of the apparent 
contradiction between the atomistic and continuum approaches to the stability of steady 
two-dimensional "simple shear". 

It seems to be an accepted result (though the exact assumptions required are unclear) that 
linear stability theory predicts stability for simple shear at all Reynolds numbers. See, for 
instance. Reference 12. On the other hand the necessary inhomogeneity due to the boundary 
conditions, and the need for thermostat and/or ergostat heat sinks, provide opportunities 
for a wide variety of much more complex models for this "simple" shear problem. 
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VII. ADDENDUM OF 24 APRIL 2012 

We wish to thank Vitaly Kuzkin for emphasizing the relevance of Denis Evans and Gary 
Morriss' paper, "Nonequilibrium Molecular Dynamics Simulation of Couette Flow in Two- 
Dimensional Fluids", Physical Review Letters 51, 1776-1779 (1983). Evans and Morriss 
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studied the shear viscosity of 896 and 3584 soft-disk particles at several strain rates and 
densities. Just as in our work in References 1-5 they observed only a small number de- 
pendence, and a nearly rate-independent viscosity rather than the divergence of viscosity 
predicted by the mode-coupling models. Their viscosity data look closely Newtonian, so 
long as the strain rate is not too large. In 1983 this was a surprise as the mode-coupling 
models of the 1970s current then predicted the divergence of viscosity in two-dimensional 
fluids. Influenced by those models, Evans and Morriss explained the observed unexpected 
stability of their own results as a "screening instability" due to transverse velocity fluctua- 
tions. By removing the three longest-wavelength fluctuations Evans and Morriss recovered 
a logarithmic rate-dependence, t] ~ ln(l/e). They commented that larger systems would 
require the removal of even more Fourier components of the transverse velocity fluctuations 
in order to squelch their observed Newtonian behavior. 
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